function [LE1, LE2] = LEs(alpha, sigma, mu, k, phi0, epsilon)
% LEs 利用 QR 分解法计算 Rulkov 双神经元及 memristor 系统前两个 Lyapunov 指数
%
% 输入:
%   alpha, sigma, mu  - Rulkov 模型参数
%   k                 - 耦合强度
%   phi0              - 忆阻器初始状态
%   epsilon           - phi 的更新系数
%
% 输出:
%   LE1 - 最大 Lyapunov 指数
%   LE2 - 第二大 Lyapunov 指数
%
% 系统状态为 [X1, Y1, X2, Y2, phi]，初始状态固定为 [0.01, 0, -0.01, 0, phi0]
% 迭代过程中，先去除瞬态 N_cut 次，再累计剩余迭代次数内 QR 分解获得的局部发散率

    % 设置迭代参数
    N_total = 20000;   % 总迭代次数
    N_cut   = 1000;    % 去除瞬态次数
    m = 5;             % 系统维数

    % 初始状态
    state = [0.01, 0, -0.01, 0, phi0];

    % 先迭代 N_cut 次去除瞬态
    for i = 1:N_cut
        state = RulkovMap(state, alpha, sigma, mu, k, epsilon);
    end

    % 初始化 QR 分解基
    Q = eye(m);
    sum_log = zeros(1, m);
    count = 0;

    % 主迭代，累计局部发散率
    for i = 1:(N_total - N_cut)
        % 计算当前状态的雅可比矩阵 (5x5)
        J = computeJacobian(state, alpha, sigma, mu, k, epsilon);
        % 传播当前正交基
        Z = J * Q;
        % 进行 QR 分解
        [Q, R] = qr(Z, 0);
        % 累计 R 对角线元素的自然对数
        for j = 1:m
            sum_log(j) = sum_log(j) + log(abs(R(j,j)));
        end
        % 更新状态
        state = RulkovMap(state, alpha, sigma, mu, k, epsilon);
        count = count + 1;
    end

    % 计算每个方向的平均发散率，即 Lyapunov 指数
    exponents = sum_log / count;
    LE1 = exponents(1);
    LE2 = exponents(2);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function next_state = RulkovMap(state, alpha, sigma, mu, k, epsilon)
% RulkovMap 根据 Rulkov 双神经元模型及 memristor 公式更新状态
%
% 输入:
%   state - 当前状态 [X1, Y1, X2, Y2, phi]
%   alpha, sigma, mu, k, epsilon - 模型参数
%
% 输出:
%   next_state - 下一时刻状态
%
% 注意：本函数中 F(x,y,alpha) 的定义假设 x<=0 区间，
% 即 F(x,y,alpha) = alpha/(1-x) + y.

    X1 = state(1); Y1 = state(2);
    X2 = state(3); Y2 = state(4);
    phi = state(5);
    
    % 对于 x<=0 区间，定义 Rulkov 映射函数
    F1 = alpha/(1 - X1) + Y1;
    F2 = alpha/(1 - X2) + Y2;
    
    % 更新规则（与 Rulkov.m 文件中一致）
    X1_next = F1 + k * (X1 - X2) * tanh(phi);
    Y1_next = Y1 - mu * (X1 + 1) + mu * sigma;
    X2_next = F2 - k * (X1 - X2) * tanh(phi);
    Y2_next = Y2 - mu * (X2 + 1) + mu * sigma;
    phi_next = phi + epsilon * (X1 - X2);
    
    next_state = [X1_next, Y1_next, X2_next, Y2_next, phi_next];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function J = computeJacobian(state, alpha, sigma, mu, k, epsilon)
% computeJacobian 根据当前状态和参数计算 Rulkov 双神经元模型的雅可比矩阵
%
% 输入:
%   state - 当前状态 [X1, Y1, X2, Y2, phi]
%   alpha, sigma, mu, k, epsilon - 模型参数
%
% 输出:
%   J - 雅可比矩阵 (5x5)
%
% 注意：此处假设系统处于 x<=0 区间，因此
% F(x,y,alpha) = alpha/(1-x)+y，
% 对应偏导数：dF/dx = alpha/(1-x)^2, dF/dy = 1.
    
    X1 = state(1); Y1 = state(2);
    X2 = state(3); Y2 = state(4);
    phi = state(5);
    
    % 对 tanh(phi) 的导数
    dtanh = sech(phi)^2;
    
    % 神经元1
    % F1 = alpha/(1-X1) + Y1
    dF1_dX1 = alpha/(1 - X1)^2;
    dF1_dY1 = 1;
    
    % X1_next = F1 + k*(X1-X2)*tanh(phi)
    dX1next_dX1 = dF1_dX1 + k * tanh(phi);
    dX1next_dY1 = dF1_dY1;
    dX1next_dX2 = - k * tanh(phi);
    dX1next_dY2 = 0;
    dX1next_dphi = k * (X1 - X2) * dtanh;
    
    % Y1_next = Y1 - mu*(X1+1) + mu*sigma
    dY1next_dX1 = - mu;
    dY1next_dY1 = 1;
    dY1next_dX2 = 0;
    dY1next_dY2 = 0;
    dY1next_dphi = 0;
    
    % 神经元2
    % F2 = alpha/(1-X2) + Y2
    dF2_dX2 = alpha/(1 - X2)^2;
    dF2_dY2 = 1;
    
    % X2_next = F2 - k*(X1-X2)*tanh(phi)
    dX2next_dX1 = - k * tanh(phi);
    dX2next_dX2 = dF2_dX2 + k * tanh(phi);
    dX2next_dY1 = 0;
    dX2next_dY2 = dF2_dY2;
    dX2next_dphi = - k * (X1 - X2) * dtanh;
    
    % Y2_next = Y2 - mu*(X2+1) + mu*sigma
    dY2next_dX1 = 0;
    dY2next_dX2 = - mu;
    dY2next_dY1 = 0;
    dY2next_dY2 = 1;
    dY2next_dphi = 0;
    
    % phi_next = phi + epsilon*(X1-X2)
    dphinext_dX1 = epsilon;
    dphinext_dX2 = - epsilon;
    dphinext_dY1 = 0;
    dphinext_dY2 = 0;
    dphinext_dphi = 1;
    
    % 组装雅可比矩阵 J (5x5)
    % 行依次对应 X1_next, Y1_next, X2_next, Y2_next, phi_next
    % 列依次对应 X1, Y1, X2, Y2, phi
    J = zeros(5,5);
    J(1,:) = [dX1next_dX1, dX1next_dY1, dX1next_dX2, dX1next_dY2, dX1next_dphi];
    J(2,:) = [dY1next_dX1, dY1next_dY1, dY1next_dX2, dY1next_dY2, dY1next_dphi];
    J(3,:) = [dX2next_dX1, dX2next_dY1, dX2next_dX2, dX2next_dY2, dX2next_dphi];
    J(4,:) = [dY2next_dX1, dY2next_dY1, dY2next_dX2, dY2next_dY2, dY2next_dphi];
    J(5,:) = [dphinext_dX1, dphinext_dY1, dphinext_dX2, dphinext_dY2, dphinext_dphi];
end
